pacman::p_load(arrow, dplyr, lubridate, ggplot2, tidyverse, sf, sfdep, sp, spNetwork, spdep, raster, spatstat, tmap, readxl, plotly)Network Kernal Density Estimation (NKDE)
1 Getting Started
osm_basemap <- tm_basemap(server = "OpenStreetMap.HOT")
imagery_basemap <- tm_basemap(server = "Esri.WorldImagery")2. Importing Data
2.1 Asaptial Data
hk_census <- read_excel("data/aspatial/hkcensus.xlsx")2.2 Geospatial Data
2.2.1 18 Districts in Hong Kong
district_18 <- st_read(dsn = "data/geospatial/hk_18Districts/",
layer = "HKDistrict18" )Reading layer `HKDistrict18' from data source
`C:\Users\Kachel Lee\ka33rina\IS415-project\Analysis\data\geospatial\hk_18Districts'
using driver `ESRI Shapefile'
Simple feature collection with 18 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 12672060 ymin: 2529946 xmax: 12739620 ymax: 2579129
Projected CRS: WGS 84 / Pseudo-Mercatorsf_district_18 <- district_18 %>% st_transform(crs = 2326)st_crs(sf_district_18)Coordinate Reference System:
User input: EPSG:2326
wkt:
PROJCRS["Hong Kong 1980 Grid System",
BASEGEOGCRS["Hong Kong 1980",
DATUM["Hong Kong 1980",
ELLIPSOID["International 1924",6378388,297,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4611]],
CONVERSION["Hong Kong 1980 Grid",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",22.3121333333333,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",114.178555555556,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",836694.05,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",819069.8,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["northing (N)",north,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["easting (E)",east,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Cadastre, engineering survey, topographic mapping (large scale)."],
AREA["China - Hong Kong - onshore and offshore."],
BBOX[22.13,113.76,22.58,114.51]],
ID["EPSG",2326]]tmap_mode("plot")
d18_map <- tm_shape(sf_district_18) +
tm_fill(col = "ENAME", title = "District", legend.show = FALSE) +
tm_borders(col = "black", lwd = 0.5) +
tm_text("ID", size = 0.5, col = "black") +
tm_layout(frame = FALSE)
d18_map
2.2.2 Recycling Points in Hong Kong
cp <- read_csv("data/aspatial/hkrecyclepoints.csv")Change the geographic location
cp_sf <- st_as_sf(cp,
coords = c("lgt","lat"),
crs = 4326) %>%
st_transform(crs= 2326)summary(cp_sf) cp_id district_id address_en address2_en
Min. : 283 Length:6520 Length:6520 Length:6520
1st Qu.: 3615 Class :character Class :character Class :character
Median : 6652 Mode :character Mode :character Mode :character
Mean : 6127
3rd Qu.: 8516
Max. :10764
waste_type legend geometry
Length:6520 Length:6520 POINT :6520
Class :character Class :character epsg:2326 : 0
Mode :character Mode :character +proj=tmer...: 0
cp_sf_1 <- cp_sf %>%
mutate(district_id = toupper(str_replace_all(district_id, "_", " ")))recycling_bins <- subset(cp_sf_1, legend == "Recycling Bins at Public Place")
recycling_spots <- subset(cp_sf_1, legend == "Recycling Spots")
private_collection_points <- subset(cp_sf_1, legend == "Private Collection Points (e.g. housing estates, shopping centres)")
ngo_collection_points <- subset(cp_sf_1, legend == "NGO Collection Points")
recycling_stations <- subset(cp_sf_1, legend == "Recycling Stations/Recycling Stores")
street_corner_recycling_shops <- subset(cp_sf_1, legend == "Street Corner Recycling Shops")
smart_bins <- subset(cp_sf_1, legend == "Smart Bin")recycling_spots_cp <- st_join(sf_district_18, recycling_spots)
ngo_cp <- st_join(sf_district_18, ngo_collection_points)
pcp_joined_data <- st_join(sf_district_18, private_collection_points)
recycling_bins_cp <- st_join(sf_district_18, recycling_bins)
recycling_stations_cp <- st_join(sf_district_18, recycling_stations)
street_corner_cp <- st_join(sf_district_18, street_corner_recycling_shops)
smart_bins_cp <- st_join(sf_district_18, smart_bins)private_collection_points_by_district <- pcp_joined_data %>%
group_by(ENAME) %>%
summarize(total_pcp = n())
ngo_cp_by_district <- ngo_cp %>%
group_by(ENAME) %>%
summarize(total_ngo_cp = n())
recycling_spots_by_district <- recycling_spots_cp %>%
group_by(ENAME) %>%
summarize(total_recycling_spots = n())
recycling_bins_by_district <- recycling_bins_cp %>%
group_by(ENAME) %>%
summarize(total_recycling_bins = n())
recycling_stations_by_district <- recycling_stations_cp %>%
group_by(ENAME) %>%
summarize(total_recycling_stations = n())
street_corner_shops_by_district <- street_corner_cp %>%
group_by(ENAME) %>%
summarize(total_street_corner = n())
smart_bins_by_district <- smart_bins_cp %>%
group_by(ENAME) %>%
summarize(total_smart_bins = n())2.3 Road Data in Hong Kong
road_data <- st_read(dsn = "data/geospatial/china-latest-free.shp",
layer = "gis_osm_roads_free_1")Transform the Data into Hong Kong Projection System
road_data_2326 <- st_transform(road_data, 2326)roads_in_hk <- st_intersection(road_data_2326, sf_district_18)write_rds(roads_in_hk, "data/rds/sf_roads_in_hk.rds")roads_in_hk <- read_rds("data/rds/sf_roads_in_hk.rds")3. Visualizing Recycling Points in Hong Kong
pcp_map <- tm_shape(private_collection_points_by_district) +
tm_fill(col = "total_pcp") +
tm_borders() +
tm_layout(legend.show = TRUE, main.title = "Distribution of Private Collection Points by District",main.title.position = "center", main.title.size = 0.75)+
tm_scale_bar()ngo_cp_map <- tm_shape(ngo_cp_by_district) +
tm_fill(col = "total_ngo_cp") +
tm_borders() +
tm_layout(legend.show = TRUE, main.title = "Distribution of NGO Points by District", main.title.position = "center", main.title.size = 0.75)+
tm_scale_bar()recycling_spots_map <- tm_shape(recycling_spots_by_district) +
tm_fill(col = "total_recycling_spots") +
tm_borders() +
tm_layout(legend.show = TRUE, main.title = "Distribution of Recycling Spots by District",main.title.position = "center", main.title.size = 0.75)+
tm_scale_bar()recycling_bins_map <- tm_shape(recycling_bins_by_district) +
tm_fill(col = "total_recycling_bins") +
tm_borders() +
tm_layout(legend.show = TRUE, main.title = "Distribution of Recycling Bins by District", main.title.position = "center", main.title.size = 0.75)+
tm_scale_bar()recycling_stations_map <- tm_shape(recycling_stations_by_district) +
tm_fill(col = "total_recycling_stations") +
tm_borders() +
tm_layout(legend.show = TRUE, main.title = "Distribution of Recycling Stations by District", main.title.position = "center", main.title.size = 0.75)+
tm_scale_bar()street_corner_shops_map <- tm_shape(street_corner_shops_by_district) +
tm_fill(col = "total_street_corner") +
tm_borders() +
tm_layout(legend.show = TRUE,main.title = "Distribution of Street Corner Shops by District", main.title.position = "center", main.title.size = 0.75)+
tm_scale_bar()smart_bin_map <- tm_shape(smart_bins_by_district) +
tm_fill(col = "total_smart_bins") +
tm_borders() +
tm_layout(legend.show = TRUE,main.title = "Distribution of Smart Bins by District",main.title.position = "center", main.title.size = 0.75)+
tm_scale_bar()tmap_arrange(pcp_map,recycling_spots_map, recycling_bins_map,ngo_cp_map, recycling_stations_map, street_corner_shops_map, smart_bin_map, ncol=2, nrow = 4)
Since we want to conduct a more accurate analysis on the recycling points, We will focus on the recycling facilities type with a higher number of data - Private Collection Points and Recycling Bins.
3.1 Private Collection Points
wm_q_pcp <- private_collection_points_by_district %>%
mutate(nb = st_contiguity(geometry, queen = TRUE),
wt = st_weights(nb,
style = "W",
allow_zero = TRUE),
.before = 1)
wm_q_pcpSimple feature collection with 18 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 801026 ymin: 801678 xmax: 863539.9 ymax: 846903.1
Projected CRS: Hong Kong 1980 Grid System
# A tibble: 18 × 5
nb wt ENAME total_pcp geometry
* <nb> <list> <chr> <int> <MULTIPOLYGON [m]>
1 <int [2]> <dbl [2]> CENTRAL & WESTERN 414 (((833047.3 816838.9, 833488…
2 <int [2]> <dbl [2]> EASTERN 200 (((843535.7 812736.3, 843530…
3 <int [1]> <dbl [1]> ISLANDS 46 (((810029 801683.4, 810019.4…
4 <int [5]> <dbl [5]> KOWLOON CITY 243 (((836281 823326.3, 836283 8…
5 <int [3]> <dbl [3]> KWAI TSING 87 (((827834 820693.1, 827834 8…
6 <int [3]> <dbl [3]> KWUN TONG 119 (((843155.9 816368.9, 843154…
7 <int [2]> <dbl [2]> NORTH 97 (((852615.5 841161.9, 852615…
8 <int [4]> <dbl [4]> SAI KUNG 163 (((840825.4 823785.4, 840827…
9 <int [7]> <dbl [7]> SHA TIN 194 (((839924.4 832031.6, 839924…
10 <int [4]> <dbl [4]> SHAM SHUI PO 190 (((834919.6 823272.4, 834921…
11 <int [3]> <dbl [3]> SOUTHERN 181 (((830176.1 815000.2, 830255…
12 <int [5]> <dbl [5]> TAI PO 121 (((846104.6 830686.6, 846265…
13 <int [6]> <dbl [6]> TSUEN WAN 92 (((821762.8 819430, 821765.7…
14 <int [2]> <dbl [2]> TUEN MUN 165 (((811619.1 831909.9, 811639…
15 <int [3]> <dbl [3]> WAN CHAI 198 (((838663.3 815002.8, 838700…
16 <int [4]> <dbl [4]> WONG TAI SIN 82 (((836530.5 823327.5, 836533…
17 <int [2]> <dbl [2]> YAU TSIM MONG 116 (((835145.2 817859.4, 835137…
18 <int [4]> <dbl [4]> YUEN LONG 178 (((811708.8 831974.4, 811719…set.seed(1234)
global_moran_perm(wm_q_pcp$total_pcp,
wm_q_pcp$nb,
wm_q_pcp$wt,
zero.policy = TRUE,
nsim = 999)
Monte-Carlo simulation of Moran I
data: x
weights: listw
number of simulations + 1: 1000
statistic = 0.16065, observed rank = 934, p-value = 0.132
alternative hypothesis: two.sidedlisa_pcp <- wm_q_pcp %>%
mutate(local_moran = local_moran(
total_pcp, nb, wt, zero.policy = TRUE, nsim = 99),
.before = 1) %>%
unnest(local_moran)lisa_sig_pcp <- lisa_pcp %>%
filter(p_ii < 0.05)
tmap_mode("plot")
tm_shape(lisa_pcp) +
tm_polygons() +
tm_borders(alpha = 0.7) +
tm_shape(lisa_pcp) +
tm_fill("mean") +
tm_borders(alpha = 0.4)
3.2 Recycling Bins
wm_q_rbins <- recycling_bins_by_district %>%
mutate(nb = st_contiguity(geometry, queen = TRUE),
wt = st_weights(nb,
style = "W",
allow_zero = TRUE),
.before = 1)
wm_q_rbinsSimple feature collection with 18 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 801026 ymin: 801678 xmax: 863539.9 ymax: 846903.1
Projected CRS: Hong Kong 1980 Grid System
# A tibble: 18 × 5
nb wt ENAME total_recycling_bins geometry
* <nb> <list> <chr> <int> <MULTIPOLYGON [m]>
1 <int [2]> <dbl [2]> CENTRAL &… 93 (((833047.3 816838.9, 83…
2 <int [2]> <dbl [2]> EASTERN 177 (((843535.7 812736.3, 84…
3 <int [1]> <dbl [1]> ISLANDS 219 (((810029 801683.4, 8100…
4 <int [5]> <dbl [5]> KOWLOON C… 92 (((836281 823326.3, 8362…
5 <int [3]> <dbl [3]> KWAI TSING 130 (((827834 820693.1, 8278…
6 <int [3]> <dbl [3]> KWUN TONG 133 (((843155.9 816368.9, 84…
7 <int [2]> <dbl [2]> NORTH 286 (((852615.5 841161.9, 85…
8 <int [4]> <dbl [4]> SAI KUNG 302 (((840825.4 823785.4, 84…
9 <int [7]> <dbl [7]> SHA TIN 214 (((839924.4 832031.6, 83…
10 <int [4]> <dbl [4]> SHAM SHUI… 110 (((834919.6 823272.4, 83…
11 <int [3]> <dbl [3]> SOUTHERN 115 (((830176.1 815000.2, 83…
12 <int [5]> <dbl [5]> TAI PO 282 (((846104.6 830686.6, 84…
13 <int [6]> <dbl [6]> TSUEN WAN 193 (((821762.8 819430, 8217…
14 <int [2]> <dbl [2]> TUEN MUN 195 (((811619.1 831909.9, 81…
15 <int [3]> <dbl [3]> WAN CHAI 121 (((838663.3 815002.8, 83…
16 <int [4]> <dbl [4]> WONG TAI … 60 (((836530.5 823327.5, 83…
17 <int [2]> <dbl [2]> YAU TSIM … 147 (((835145.2 817859.4, 83…
18 <int [4]> <dbl [4]> YUEN LONG 327 (((811708.8 831974.4, 81…set.seed(1234)
global_moran_perm(wm_q_rbins$total_recycling_bins,
wm_q_rbins$nb,
wm_q_rbins$wt,
zero.policy = TRUE,
nsim = 999)
Monte-Carlo simulation of Moran I
data: x
weights: listw
number of simulations + 1: 1000
statistic = 0.48951, observed rank = 998, p-value = 0.004
alternative hypothesis: two.sidedlisa_rbins <- wm_q_rbins %>%
mutate(local_moran = local_moran(
total_recycling_bins, nb, wt, zero.policy = TRUE, nsim = 99),
.before = 1) %>%
unnest(local_moran)lisa_sig_rbins <- lisa_rbins %>%
filter(p_ii < 0.05)
tmap_mode("plot")
tm_shape(lisa_rbins) +
tm_polygons() +
tm_borders(alpha = 0.7) +
tm_shape(lisa_rbins) +
tm_fill("mean") +
tm_borders(alpha = 0.4)
Considering the Hot-Spot Analysis on the Private Collection Points and Public Recycling bins, they are maininly locate far from the residential area. Therefore, we decide to refer to the census data in Hong Kong.
ggplot(hk_census, aes(x = DNAME, y = POPULATION)) +
geom_bar(stat = "identity", fill = "steelblue", color = "black") +
geom_text(aes(label = POPULATION), vjust = -0.5) +
theme_minimal() +
labs(title = "Population by District",
x = "District Name",
y = "Population") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
Since Sha Tin has the most population, we will focus on Sha Tin district. Moreover, we want to choose another district with both commercial and residential area. Therefore, we will be analyzing on Yau Tsim Mong Area as well.
sha_tin_recycling_points <- cp_sf_1 %>%
filter(district_id == "SHA TIN" &
(legend == "Private Collection Points (e.g. housing estates, shopping centres)" | legend == "Recycling Bins at Public Place"))ytm_recycling_points <- cp_sf_1 %>%
filter(district_id == "YAU TSIM MONG" &
(legend == "Private Collection Points (e.g. housing estates, shopping centres)" | legend == "Recycling Bins at Public Place"))wanchai_recycling_points <- cp_sf_1 %>%
filter(district_id == "WAN CHAI" &
(legend == "Private Collection Points (e.g. housing estates, shopping centres)" | legend == "Recycling Bins at Public Place"))4. Sha Tin District
shatin <- sf_district_18 %>%
filter(ENAME == "SHA TIN")plot(st_geometry((shatin)))
roads_in_shatin <- st_intersection(roads_in_hk, shatin)write_rds(roads_in_shatin, "data/rds/roads_in_shatin.rds")roads_in_shatin <- read_rds("data/rds/roads_in_shatin.rds")private_cp_shatin <- st_intersection(private_collection_points, shatin)public_cp_shatin <- st_intersection(recycling_bins, shatin)cp_shatin <- st_intersection(sha_tin_recycling_points, shatin)write_rds(private_cp_shatin, "data/rds/private_cp_shatin.rds")write_rds(public_cp_shatin, "data/rds/public_cp_shatin.rds")write_rds(cp_shatin, "data/rds/cp_shatin.rds")private_cp_shatin <- read_rds("data/rds/private_cp_shatin.rds")public_cp_shatin <- read_rds("data/rds/public_cp_shatin.rds")cp_shatin <- read_rds("data/rds/cp_shatin.rds")plot(st_geometry(private_cp_shatin))
plot(st_geometry(public_cp_shatin))
plot(st_geometry(cp_shatin))
unique_types <- unique(st_geometry_type(roads_in_shatin))
unique_types[1] POINT GEOMETRYCOLLECTION LINESTRING MULTILINESTRING
[5] MULTIPOINT
18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLEroads_in_shatin <- roads_in_shatin %>%
filter(st_geometry_type(geometry) %in% c("LINESTRING", "MULTILINESTRING"))if ("LINESTRING" %in% unique_types) {
roads_in_shatin <- st_cast(roads_in_shatin, "LINESTRING")
} else {
# handle the case when no linestrings are found
stop("No linestrings found in roads_in_shatin")
}tmap_mode('plot')
tm_shape(roads_in_shatin, geometry_type = "lines") +
tm_lines()
tmap_mode('plot')
# Plot roads
tm_shape(roads_in_shatin, geometry_type = "lines") +
tm_lines(lwd = 1, col = "blue") +
# Plot Punggol area boundary
tm_shape(shatin) +
tm_borders() +
# Plot Grab origin data
tm_shape(private_cp_shatin) +
tm_dots()
tmap_mode("view")
# Create the map using tmap
cp_shatin_map <- osm_basemap +
tm_shape(cp_shatin) +
tm_dots(col = "legend", palette = c("blue", "red"), border.col = "black", size = 0.05) +
tm_layout(title = "Recycling Bin Locations in Sha Tin")
cp_shatin_maproads_lines_shatin <- roads_in_shatin[st_geometry_type(roads_in_shatin) == "LINESTRING", ]write_rds(roads_lines_shatin, "data/rds/roads_lines_shatin.rds")roads_lines_shatin <- read_rds("data/rds/roads_lines_shatin.rds")# Apply lixelize_lines with mindist
lixels_shatin <- lixelize_lines(roads_lines_shatin,5000, mindist = 2500)samples_shatin <- lines_center(roads_lines_shatin)4.1 Private Collection Points in Sha Tin
densities_shatin_private <- nkde(roads_lines_shatin,
events = private_cp_shatin,
w = rep(1,nrow(private_cp_shatin)),
samples = samples_shatin,
kernel_name = "quartic",
bw = 300,
div= "bw",
method = "simple",
digits = 1,
tol = 1,
grid_shape = c(1,1),
max_depth = 8,
agg = 5, #we aggregate events within a 5m radius (faster calculation)
sparse = TRUE,
verbose = FALSE)write_rds(densities_shatin_private , "data/rds/densities_shatin_private.rds")densities_shatin_private <- read_rds("data/rds/densities_shatin_private.rds")samples_shatin$density_private <- densities_shatin_private
lixels_shatin$density_private <- densities_shatin_private# rescaling to help the mapping
samples_shatin$density_private <- samples_shatin$density_private*1000
lixels_shatin$density_private <- lixels_shatin$density_private*1000tmap_mode('view')
shatin_density_private <- osm_basemap+
tm_shape(lixels_shatin)+
tm_lines(col="density_private")+
tm_shape(private_cp_shatin)+
tm_dots()
shatin_density_private4.2 Public Recycling Bins in Sha Tin
densities_shatin_public <- nkde(roads_lines_shatin,
events = public_cp_shatin,
w = rep(1,nrow(public_cp_shatin)),
samples = samples_shatin,
kernel_name = "quartic",
bw = 300,
div= "bw",
method = "simple",
digits = 1,
tol = 1,
grid_shape = c(1,1),
max_depth = 8,
agg = 5, #we aggregate events within a 5m radius (faster calculation)
sparse = TRUE,
verbose = FALSE)write_rds(densities_shatin_public , "data/rds/densities_shatin_public.rds")densities_shatin_public <- read_rds("data/rds/densities_shatin_public.rds")samples_shatin$density_public <- densities_shatin_public
lixels_shatin$density_public <- densities_shatin_public# rescaling to help the mapping
samples_shatin$density_public <- samples_shatin$density_public*1000
lixels_shatin$density_public <- lixels_shatin$density_public*1000tmap_mode('view')
shatin_density_public <- osm_basemap+
tm_shape(lixels_shatin)+
tm_lines(col="density_public")+
tm_shape(public_cp_shatin)+
tm_dots()
shatin_density_public4.3 Both Private and Public Bins in Sha Tin Districts
densities_shatin_ppcp <- nkde(roads_lines_shatin,
events = cp_shatin,
w = rep(1,nrow(cp_shatin)),
samples = samples_shatin,
kernel_name = "quartic",
bw = 300,
div= "bw",
method = "simple",
digits = 1,
tol = 1,
grid_shape = c(1,1),
max_depth = 8,
agg = 5, #we aggregate events within a 5m radius (faster calculation)
sparse = TRUE,
verbose = FALSE)write_rds(densities_shatin_ppcp , "data/rds/densities_shatin_ppcp.rds")densities_shatin_ppcp <- read_rds("data/rds/densities_shatin_ppcp.rds")samples_shatin$density_ppcp <- densities_shatin_ppcp
lixels_shatin$density_ppcp <- densities_shatin_ppcp# rescaling to help the mapping
samples_shatin$density_ppcp <- samples_shatin$density_ppcp*1000
lixels_shatin$density_ppcp <- lixels_shatin$density_ppcp*1000tmap_mode('view')
shatin_density_ppcp <- osm_basemap +
tm_shape(lixels_shatin) +
tm_lines(col = "density_ppcp") +
tm_shape(cp_shatin) +
tm_layout(
legend.position = c("left", "top"),
legend.text.size = 0.5, # Adjust this value as needed for smaller legend text
legend.title.size = 0.6 # Adjust this value as needed for smaller legend title
)+
tm_dots(col = "legend", palette = c("blue", "red"), border.col = "black", size = 0.05)
shatin_density_ppcpkfun_shatin_ppcp <- kfunctions(roads_in_shatin,
cp_shatin,
start = 0,
end = 1000,
step = 50,
width = 50,
nsim = 50,
resolution = 50,
verbose = FALSE,
conf_int = 0.05,
agg = 25
)write_rds(kfun_shatin_ppcp, "data/rds/kfun_shatin_ppcp.rds")kfun_shatin_ppcp <- read_rds("data/rds/kfun_shatin_ppcp.rds")kfun_shatin_ppcp$plotk
5. Yau Tsim Mong Area
ytm <- sf_district_18 %>%
filter(ENAME == "YAU TSIM MONG")plot(st_geometry((ytm)))
roads_in_ytm <- st_intersection(roads_in_hk, ytm)write_rds(roads_in_ytm, "data/rds/roads_in_ytm.rds")roads_in_ytm <- read_rds("data/rds/roads_in_ytm.rds")private_cp_ytm <- st_intersection(private_collection_points, ytm)public_cp_ytm <- st_intersection(recycling_bins, ytm)cp_ytm <- st_intersection(ytm_recycling_points, ytm)write_rds(private_cp_ytm, "data/rds/private_cp_ytm.rds")write_rds(public_cp_ytm, "data/rds/public_cp_ytm.rds")write_rds(cp_ytm, "data/rds/cp_ytm.rds")private_cp_ytm <- read_rds("data/rds/private_cp_ytm.rds")public_cp_ytm <- read_rds("data/rds/public_cp_ytm.rds")cp_ytm <- read_rds("data/rds/cp_ytm.rds")plot(st_geometry(private_cp_ytm))
plot(st_geometry(public_cp_ytm))
unique_types <- unique(st_geometry_type(roads_in_ytm))
unique_types[1] LINESTRING MULTILINESTRING POINT GEOMETRYCOLLECTION
18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLEroads_in_ytm <- roads_in_ytm %>%
filter(st_geometry_type(geometry) %in% c("LINESTRING", "MULTILINESTRING"))if ("LINESTRING" %in% unique_types) {
roads_in_ytm <- st_cast(roads_in_ytm, "LINESTRING")
} else {
# handle the case when no linestrings are found
stop("No linestrings found in roads_in_ytm")
}tmap_mode('plot')
tm_shape(roads_in_ytm, geometry_type = "lines") +
tm_lines()5.1 Private Collection Poitns in YTM
tmap_mode('plot')
# Plot roads
tm_shape(roads_in_ytm, geometry_type = "lines") +
tm_lines(lwd = 1, col = "blue") +
tm_shape(ytm) +
tm_borders() +
tm_shape(private_cp_ytm) +
tm_dots()tmap_mode("view")
# Create the map using tmap
cp_ytm_map <- osm_basemap +
tm_shape(cp_ytm) +
tm_dots(col = "legend", palette = c("blue", "red"), border.col = "black", size = 0.05) +
tm_layout(title = "Recycling Bin Locations in Yau Tsim Mong")
cp_shatin_maproads_lines_ytm <- roads_in_ytm[st_geometry_type(roads_in_ytm) == "LINESTRING", ]write_rds(roads_lines_ytm, "data/rds/roads_lines_ytm.rds")roads_lines_ytm <- read_rds("data/rds/roads_lines_ytm.rds")# Apply lixelize_lines with mindist
lixels_ytm <- lixelize_lines(roads_lines_ytm,5000, mindist = 2500)samples_ytm <- lines_center(lixels_ytm)densities_ytm_private <- nkde(roads_lines_ytm,
events = private_cp_ytm,
w = rep(1,nrow(private_cp_ytm)),
samples = samples_ytm,
kernel_name = "quartic",
bw = 300,
div= "bw",
method = "simple",
digits = 1,
tol = 1,
grid_shape = c(1,1),
max_depth = 8,
agg = 5, #we aggregate events within a 5m radius (faster calculation)
sparse = TRUE,
verbose = FALSE)write_rds(densities_ytm_private , "data/rds/densities_ytm_private.rds")densities_ytm_private <- read_rds("data/rds/densities_ytm_private.rds")samples_ytm$density_private <- densities_ytm_private
lixels_ytm$density_private <- densities_ytm_private# rescaling to help the mapping
samples_ytm$density_private <- samples_ytm$density_private*1000
lixels_ytm$density_private <- lixels_ytm$density_private*1000tmap_mode('view')
ytm_density_private <- osm_basemap +
tm_shape(lixels_ytm)+
tm_lines(col="density_private")+
tm_shape(private_cp_ytm)+
tm_dots()
ytm_density_private5.2 Public Recycling Bins in YTM
densities_ytm_public <- nkde(roads_lines_ytm,
events = public_cp_ytm,
w = rep(1,nrow(public_cp_ytm)),
samples = samples_ytm,
kernel_name = "quartic",
bw = 300,
div= "bw",
method = "simple",
digits = 1,
tol = 1,
grid_shape = c(1,1),
max_depth = 8,
agg = 5, #we aggregate events within a 5m radius (faster calculation)
sparse = TRUE,
verbose = FALSE)write_rds(densities_ytm_public , "data/rds/densities_ytm_public.rds")densities_ytm_public <- read_rds("data/rds/densities_ytm_public.rds")samples_ytm$density_public <- densities_ytm_public
lixels_ytm$density_public <- densities_ytm_public# rescaling to help the mapping
samples_ytm$density_public <- samples_ytm$density_public*1000
lixels_ytm$density_public <- lixels_ytm$density_public*1000tmap_mode('view')
ytm_density_public <-osm_basemap +
tm_shape(lixels_ytm)+
tm_lines(col="density_public")+
tm_shape(private_cp_ytm)+
tm_dots()
ytm_density_public5.3 Both Public and Private Recycling Bins in YTM
densities_ytm_ppcp <- nkde(roads_lines_ytm,
events = cp_ytm,
w = rep(1,nrow(cp_ytm)),
samples = samples_ytm,
kernel_name = "quartic",
bw = 300,
div= "bw",
method = "simple",
digits = 1,
tol = 1,
grid_shape = c(1,1),
max_depth = 8,
agg = 5, #we aggregate events within a 5m radius (faster calculation)
sparse = TRUE,
verbose = FALSE)write_rds(densities_ytm_ppcp , "data/rds/densities_ytm_ppcp.rds")densities_ytm_ppcp <- read_rds("data/rds/densities_ytm_ppcp.rds")samples_ytm$density_ppcp <- densities_ytm_ppcp
lixels_ytm$density_ppcp <- densities_ytm_ppcp# rescaling to help the mapping
samples_ytm$density_ppcp <- samples_ytm$density_ppcp*1000
lixels_ytm$density_ppcp <- lixels_ytm$density_ppcp*1000tmap_mode('view')
ytm_density_ppcp <-osm_basemap +
tm_shape(lixels_ytm)+
tm_lines(col="density_ppcp")+
tm_shape(cp_ytm)+
tm_dots(col = "legend", palette = c("blue", "red"), border.col = "black", size = 0.05)
ytm_density_ppcpkfun_ytm_ppcp <- kfunctions(roads_in_ytm,
cp_ytm,
start = 0,
end = 1000,
step = 50,
width = 50,
nsim = 50,
resolution = 50,
verbose = FALSE,
conf_int = 0.05,
agg = 20
)write_rds(kfun_ytm_ppcp, "data/rds/kfun_ytm_ppcp.rds")kfun_ytm_ppcp <- read_rds("data/rds/kfun_ytm_ppcp.rds")kfun_ytm_ppcp$plotk
6. Wan Chai (Hong Kong Island)
wanchai <- sf_district_18 %>%
filter(ENAME == "WAN CHAI")plot(st_geometry((wanchai)))
roads_in_wanchai <- st_intersection(roads_in_hk, wanchai)write_rds(roads_in_wanchai, "data/rds/roads_in_wanchai.rds")roads_in_wanchai <- read_rds("data/rds/roads_in_wanchai.rds")private_cp_wanchai <- st_intersection(private_collection_points, wanchai)public_cp_wanchai<- st_intersection(recycling_bins, wanchai)cp_wanchai <- st_intersection(wanchai_recycling_points, wanchai)write_rds(private_cp_wanchai, "data/rds/private_cp_wanchai.rds")write_rds(public_cp_wanchai, "data/rds/public_cp_wanchai.rds")write_rds(cp_wanchai, "data/rds/cp_wanchai.rds")private_cp_wanchai <- read_rds("data/rds/private_cp_wanchai.rds")public_cp_wanchai <- read_rds("data/rds/public_cp_wanchai.rds")cp_wanchai <- read_rds("data/rds/cp_wanchai.rds")plot(st_geometry(private_cp_wanchai))
plot(st_geometry(cp_wanchai))
unique_types <- unique(st_geometry_type(roads_in_wanchai))
unique_types[1] GEOMETRYCOLLECTION POINT LINESTRING MULTIPOINT
[5] MULTILINESTRING
18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLEroads_in_wanchai <- roads_in_wanchai %>%
filter(st_geometry_type(geometry) %in% c("LINESTRING", "MULTILINESTRING"))if ("LINESTRING" %in% unique_types) {
roads_in_wanchai <- st_cast(roads_in_wanchai, "LINESTRING")
} else {
# handle the case when no linestrings are found
stop("No linestrings found in roads_in_wanchai")
}tmap_mode('plot')
tm_shape(roads_in_wanchai, geometry_type = "lines") +
tm_lines()6.3 Both Private and Public Bins in Wan Chai Districts
roads_lines_wanchai <- roads_in_wanchai[st_geometry_type(roads_in_wanchai) == "LINESTRING", ]write_rds(roads_lines_wanchai, "data/rds/roads_lines_wanchai.rds")roads_lines_wanchai <- read_rds("data/rds/roads_lines_wanchai.rds")# Apply lixelize_lines with mindist
lixels_wanchai <- lixelize_lines(roads_lines_wanchai,5000, mindist = 2500)samples_wanchai <- lines_center(lixels_wanchai)densities_wanchai_ppcp <- nkde(roads_lines_wanchai,
events = cp_wanchai,
w = rep(1,nrow(cp_wanchai)),
samples = samples_wanchai,
kernel_name = "quartic",
bw = 300,
div= "bw",
method = "simple",
digits = 1,
tol = 1,
grid_shape = c(1,1),
max_depth = 8,
agg = 5, #we aggregate events within a 5m radius (faster calculation)
sparse = TRUE,
verbose = FALSE)write_rds(densities_wanchai_ppcp , "data/rds/densities_wanchai_ppcp.rds")densities_wanchai_ppcp <- read_rds("data/rds/densities_wanchai_ppcp.rds")samples_wanchai$density_ppcp <- densities_wanchai_ppcp
lixels_wanchai$density_ppcp <- densities_wanchai_ppcp# rescaling to help the mapping
samples_wanchai$density_ppcp <- samples_wanchai$density_ppcp*1000
lixels_wanchai$density_ppcp <- lixels_wanchai$density_ppcp*1000tmap_mode('view')
wanchai_density_ppcp <- osm_basemap +
tm_shape(lixels_wanchai) +
tm_lines(col = "density_ppcp") +
tm_shape(cp_wanchai) +
tm_layout(
legend.position = c("left", "top"),
legend.text.size = 0.5, # Adjust this value as needed for smaller legend text
legend.title.size = 0.6 # Adjust this value as needed for smaller legend title
)+
tm_dots(col = "legend", palette = c("blue", "red"), border.col = "black", size = 0.05)
wanchai_density_ppcpkfun_wanchai_ppcp <- kfunctions(roads_in_wanchai,
cp_wanchai,
start = 0,
end = 1000,
step = 50,
width = 50,
nsim = 50,
resolution = 50,
verbose = FALSE,
conf_int = 0.05,
agg = 15
)write_rds(kfun_wanchai_ppcp, "data/rds/kfun_wanchai_ppcp.rds")kfun_wanchai_ppcp <- read_rds("data/rds/kfun_wanchai_ppcp.rds")kfun_wanchai_ppcp$plotk